% Obtain equilibrium domestic and export shares and other moments from
% distribution of efficiency and re-scaled fixed costs, as well as grid-level optimal domestic share matricies 
function [sD,sX,EqMoments,omega,TotalExports,TotalImports] = SolveModel(phi,fc,fx,FXM,FX,FM,SDM, SDXM, zvec, yvec, param)

betastuff = param.betastuff ; betastuff2 =  param.betastuff2 ; betastuff3 = param.betastuff3 ;
eps = param.eps;
g = param.gamma; 
sg = param.sg;
eta = param.eta;
theta = param.theta;
beta = param.beta;

zz = fc./phi.^(param.sg-1); 
yy = fx./phi.^(param.sg-1); 

%Cap export costs;
yy = min(yy,max(yvec));
yy = max(yy,min(yvec));

if max(zz) > max(zvec)
    warning('Interpolation of z out of range')
end
if max(yy) > max(yvec)
    warning('Interpolation of y out of range')
end

% Interpolate to population of firms characterized by phi, fc and fx.

sdM = interp1(zvec,SDM,zz,'spline');
sdM = min(1-1e-08,sdM);
sdM = max(1e-08,sdM);

sdXM=interp2(zvec,yvec,SDXM,zz,yy,'spline');
sdXM = min(1-1e-08,sdXM);
sdXM = max(1e-08,sdXM);

% International status
% Compute re-scaled profits for each status

piD = betastuff *  phi.^(sg-1); % purely domestic (non-importer and non-exporter)

piM2 = betastuff * phi.^(sg-1).* sdM.^(g *(1-sg)/(eps-1))  - ...
        betastuff3* (sdM.^(-1)-1).^(1/(eta*(eps-1))) .* fc * eta*g*(sg-1); % importer-only

piX2 = betastuff *  phi.^(sg-1) + betastuff2 * phi.^(theta*(sg-1)) .* fx.^(1-theta) / theta ; % exporter-only
    
piXM2 = betastuff *  phi.^(sg-1).* sdXM.^(-g *(sg-1)/(eps-1)) + ...
       betastuff2 * phi.^(theta*(sg-1)) .* sdXM.^(-theta* g *(sg-1)/(eps-1)) .* fx.^(1-theta) / theta  - ...
      (sdXM.^(-1)-1).^(1/(eta*(eps-1))) .* fc * eta*g*(sg-1)*betastuff3 ; % global
     
PI = [piD piM2-FM piX2-FX piXM2-FXM]; % include fixed costs to overall international strategy

checkPI = isnan(PI);

if max(max(checkPI))==1
    warning('NAN entries in PI')
end
if isreal(PI)==0
    warning('Complex entries in PI')        
end

[~,In] = max(PI,[],2); % In gives optimal international status with, 
% 1 denoting purely-domestic, 2 denoting importer-only, 3 denoting 
% exporter-only, and 4 denoting global


S = In;
    
% Optimal domestic shares:
sD = (S==1) + (S==2) .* sdM + (S==3) + (S==4) .*sdXM;

% Associated export share:
sX = (1 +  phi.^(-(sg-1)*(theta-1)) .* sD .^(g / (eps-1)* (sg-1)*(theta-1)) .*  beta ^ (-g * eps / (eps-1) * (sg-1) * (theta-1)) .* fx.^(theta-1)).^(-1) ;
sX = ((S==3) + (S==4)) .* sX;

% Moments
% Variable Profits (i.e., profit excluding any fixed costs):
vaD = piD; % purely domestic
vaM = betastuff * phi.^(sg-1).* sdM.^(g *(1-sg)/(eps-1)); % importer-only
vaX =  betastuff *  phi.^(sg-1) + ...
       betastuff2 * phi.^(theta*(sg-1)) .* fx.^(1-theta); % expoter-only
vaXM=  betastuff *  phi.^(sg-1).* sdXM.^(g *(1-sg)/(eps-1)) + ...
       betastuff2 * phi.^(theta*(sg-1)) .* sdXM.^(theta* g *(1-sg)/(eps-1)) .* fx.^(1-theta); % global

va = (S==1) .* vaD + (S==2) .* vaM + (S==3) .* vaX + (S==4) .* vaXM;
% Material shares:
omega = va/sum(va); % materials and sales are both proportional to variable profits

% Fractions of importers, exporters and importer-exporters
Firms_Importers = sum ((S==2))/length(sD);
Firms_Exporters = sum ((S==3))/length(sD);
Firms_ImpExp    = sum ((S==4))/length(sD);
 
% Aggregate import and export Share
Aggsi = sum( (1-sD) .* omega );
AggsX = sum( sX .* omega );

% Standard deviations of log materials, import shares, and export shares:
STD_lnva = std(log(va));
STD_si = (var(1-sD))^(1/2);
STD_xi = (var(sX))^(1/2);

% Pairwise correlations between log materials, import shares, and export shares:
corr_sixi = corr(1-sD,sX);
corr_lnvasi = corr(log(va),1-sD);
corr_lnvaxi = corr(log(va),sX);


% Exports (re-scaled) in domestic labor:  
xX = sg* betastuff2 * phi.^(theta*(sg-1)) .* fx.^(1-theta); % exporter-only
xXM = sg* betastuff2 * phi.^(theta*(sg-1)) .* sdXM.^(theta* g *(1-sg)/(eps-1)) .* fx.^(1-theta); % global
TotalExports = mean( (S==3) .* xX + (S==4) .* xXM );

% Imports of materials (use that materials are proportional to variable
% profits):
mM = vaM * g .* (1-sdM) * (sg-1); % importer-only
mXM= vaXM * g .* (1-sdXM) * (sg-1); % global
TotalImports = mean( (S==2) .* mM + (S==4) .* mXM );

% Save results

EqMoments.Aggsi = Aggsi;
EqMoments.AggsX = AggsX;
EqMoments.STD_lnva = STD_lnva;
EqMoments.STD_si = STD_si;
EqMoments.STD_xi = STD_xi;
EqMoments.corr_sixi = corr_sixi;
EqMoments.corr_lnvasi = corr_lnvasi;
EqMoments.corr_lnvaxi = corr_lnvaxi;
EqMoments.Firms_Importers = Firms_Importers;
EqMoments.Firms_Exporters = Firms_Exporters;
EqMoments.Firms_ImpExp = Firms_ImpExp;

end